{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing frame dragging in Kerr space-time\n",
    "\n",
    "### Importing required modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "import numpy as np\n",
    "from einsteinpy.metric import Kerr\n",
    "from einsteinpy.coordinates import BoyerLindquistDifferential\n",
    "from einsteinpy.bodies import Body\n",
    "from einsteinpy.geodesic import Geodesic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining position/velocity of test particle\n",
    "  - Initial velocity is kept 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 1.989e30 * u.kg\n",
    "a = 0.3 * u.m\n",
    "BL_obj = BoyerLindquistDifferential(49.95e5 * u.km, np.pi / 2 * u.rad, np.pi * u.rad,\n",
    "                                    0 * u.km / u.s, 0 * u.rad / u.s, 0 * u.rad / u.s,\n",
    "                                    a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "end_lambda = ((1 * u.year).to(u.s)).value / 930\n",
    "# Choosing stepsize for ODE solver to be 0.02 minutes\n",
    "stepsize = ((0.02 * u.min).to(u.s)).value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "obj = Kerr.from_coords(BL_obj, M)\n",
    "ans = obj.calculate_trajectory(\n",
    "    end_lambda=end_lambda, OdeMethodKwargs={\"stepsize\": stepsize}, return_cartesian=True\n",
    ")[1]\n",
    "x, y = ans[:,1], ans[:,2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting the trajectory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEJCAYAAAC+I6F6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGnVJREFUeJzt3XuQ1eWd5/H3p7sBRZRrc5GLDYgiONGNZ1BjbUpFFLOjqBNrMWZldpywpuJOprYyG11qo6tFRWd2151czAyTMVEHNa5VE0iMEkSNqR0RmkRU0JYGBFoQW8ArCnb3d/84P8yhOU03/Zzu05fPq+rU+f2e8/ye3/M0l0//nt/lKCIwMzPrrIpyd8DMzHo3B4mZmSVxkJiZWRIHiZmZJXGQmJlZEgeJmZklcZCYmVkSB4mZmSVxkJiZWZKqcnegO4waNSpqamrK3Q0zs15l3bp170REdXv1+kWQ1NTUUFtbW+5umJn1KpK2daSep7bMzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySOEjMzCxJSYJE0lxJdZLqJd1S5PNBkn6Wff6CpJqCz27NyuskXVZQfp+ktyW90qqtEZJWStqUvQ8vxRjMzKxzkoNEUiXwQ+ByYAZwnaQZrardCOyLiFOBe4C7s21nAPOBmcBc4N6sPYCfZmWt3QKsiohpwKps3czMyqQURySzgPqI2BIRB4FHgHmt6swD7s+WHwNmS1JW/khEHIiIrUB91h4R8Rywt8j+Ctu6H7iqBGMwM7NOKkWQjAd2FKw3ZGVF60REE/AeMLKD27Y2JiJ2ZW3tAkYXqyRpoaRaSbWNjY0dHIqZmR2rUgSJipRFB+t0ZNtOiYglEZGLiFx1dbvPHDMzs04qRZA0ABML1icAO9uqI6kKGEp+2qoj27a2W9K4rK1xwNud7rmZmSUrRZCsBaZJmixpIPmT58tb1VkOLMiWvww8HRGRlc/PruqaDEwD1rSzv8K2FgDLSjAGMzPrpOQgyc553AysAF4FHo2IDZLukHRlVu2fgJGS6oH/QnalVURsAB4FNgJPAt+IiGYASQ8DzwOnS2qQdGPW1l3AHEmbgDnZupmZlYnyBwZ9Wy6XC38fiZnZsZG0LiJy7dXzne1mZpbEQWJmZkkcJGZmlsRBYmZmSRwkZmaWxEFiZmZJHCRmZpbEQWJmZkkcJGZmlsRBYmZmSRwkZmaWxEFiZmZJHCRmZpbEQWJmZkkcJGZmlsRBYmZmSRwkZmaWxEFiZmZJHCRmZpbEQWJmZkkcJGZmlsRBYmZmSRwkZmaWxEFiZmZJHCRmZpbEQWJmZklKEiSS5kqqk1Qv6ZYinw+S9LPs8xck1RR8dmtWXifpsvbalPRTSVslvZi9zi7FGMzMrHOqUhuQVAn8EJgDNABrJS2PiI0F1W4E9kXEqZLmA3cD/17SDGA+MBM4GXhK0mnZNkdr868j4rHUvpuZWbpSHJHMAuojYktEHAQeAea1qjMPuD9bfgyYLUlZ+SMRcSAitgL1WXsdadPMzHqAUgTJeGBHwXpDVla0TkQ0Ae8BI4+ybXttLpb0kqR7JA0qwRjMzKyTShEkKlIWHaxzrOUAtwLTgT8GRgDfLtopaaGkWkm1jY2NxaqYmVkJlCJIGoCJBesTgJ1t1ZFUBQwF9h5l2zbbjIhdkXcA+An5abAjRMSSiMhFRK66urqTQzMzs/aUIkjWAtMkTZY0kPzJ8+Wt6iwHFmTLXwaejojIyudnV3VNBqYBa47WpqRx2buAq4BXSjAGMzPrpOSrtiKiSdLNwAqgErgvIjZIugOojYjlwD8BD0qqJ38kMj/bdoOkR4GNQBPwjYhoBijWZrbLpZKqyU9/vQjclDoGMzPrPOUPDPq2XC4XtbW15e6GmVmvImldROTaq+c7283MLImDxMzMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLImDxMzMkpQkSCTNlVQnqV7SLUU+HyTpZ9nnL0iqKfjs1qy8TtJl7bUpaXLWxqaszYGlGIOZmXVOcpBIqgR+CFwOzACukzSjVbUbgX0RcSpwD3B3tu0MYD4wE5gL3Cupsp027wbuiYhpwL6sbTOzfm/p0qXU1NRQUVFBTU0NS5cu7Zb9luKIZBZQHxFbIuIg8Agwr1WdecD92fJjwGxJysofiYgDEbEVqM/aK9pmts3FWRtkbV5VgjGYmfVqS5cuZeHChWzbto2IYNu2bSxcuLBbwqSqBG2MB3YUrDcA57ZVJyKaJL0HjMzKV7fadny2XKzNkcC7EdFUpH6XqLnl8a5svserAFqAyoIyZevNBe9VRT5vyparBBJEQFPAgAoYUAmfNmf1BVUV0NSSLz8a6Q/bDqzKt3loua36A6vEgU/js+WDTcGgAfnfoSLisPU/bCeOq6rg40+bj1g+fkAlEcEnTS0cn3W49bokBg+o4KODLUAcVrb/02DwAB32fsLACiJg/6ctnDCwIvvJHe3nIE4YWMmHB5o+W/7oYPMR70MGVRERny0fTWVlJedNrear59dQWdnOH4T1OIsWLWL//v2Hle3fv59FixZx/fXXd+m+SxEkxf7GRwfrtFVe7EjpaPWP7JS0EFgIMGnSpGJVrANasvfmVuVNbbwXFRz2p/RxS0HDhzS3eu+IA20sF+1A6+XWHWi93rozHVkutt67PPHqHh564Q3OnjCU0UNP4JuXTKOqqhT/TVhX2759+zGVl1Ip/oY0ABML1icAO9uo0yCpChgK7G1n22Ll7wDDJFVlRyXF9gVARCwBlgDkcrmiYWPt8xFJ/zkiefO9g2zYvZ+6xk+oa/wEgF+8uINZNcMZMGAgt185gwEDBhy1DSufSZMmsW3btqLlXa0UQbIWmCZpMvAm+ZPnX2lVZzmwAHge+DLwdESEpOXAQ5L+N3AyMA1YQ/5f0RFtZts8k7XxSNbmshKMoU1v3PXvurJ5sx6jubmZB/7fFl7YuoePDjTz2y3v8sa7n/LGi28D8Myrb3HB1OF8eDA4b+ooT4H1MIsXL2bhwoWHTW8NHjyYxYsXd/m+FZH+y7qkLwH/h/wvovdFxGJJdwC1EbFc0nHAg8C/IX8kMj8itmTbLgL+nPwvsH8VEU+01WZWPoV8iIwAfg98NSKOOrGRy+WitrY2eZxm/UVLSwu/fmUny37/JhEtPFW3l6ZW/1WcNW4wY086jmGDB3garIdYunQpixYtYvv27UyaNInFixcnnR+RtC4icu3WK0WQ9HQOErPOiwhefXMfP3h6EycMEFv2fMK6hg+OqFczbCCzaoYBorKykovPGMPsGWOpqPB9z72Vg6SAg8SsdFpaWli5YRdPb3yLd/cfpKJCPFW3h0+LXK9wQc1Qxg877rP1qqoqbvhCDaePG0r+an7ryRwkBRwkZl2n8IhlSHahwNa9H7N2+/tF6wu4bPpIakYOZu/+gwhRfdLxXHH2yZw21gHTkzhICjhIzLpXc3MzD/7rVl7d9R7R8of/Yzbu/ohXdn1YdJtKYM70kUypHsLY4YO5/txTfDK/zBwkBRwkZj1DU1MT9/y6jnc+/IThxw/47Ihky96PWbvtvcPqXlAzlAkjjmfkCYMcLGXiICngIDHr2QoD5sMDzTy+ofGIOp+fcCLn1gwH5c+1fHO2rxLrag6SAg4Ss96jpaWFpza+xYvb9rLnowNUqIJ/Wb+LT1pdf/z58UOYNmYIo4Ycx9mnDGf2Gb5CrNQcJAUcJGa9V0Tw2s53efD5Nxh6XBUIare/z9pt7x5W7/IZ1dSMPN5HKyXkICngIDHrW9qbCvvi1OGcOf4kAMYO8/mVzupokDiyzazXqaqq4q+/NBPIT4XN2/gW67fvoyVaePWtj3h20x6e27zvs/pPbdjNRdNHUVVVxVdmTXKolJiPSMysT2lubuafn3+D3e9/zL79n7Js/VvsL7hb8otTRzhUOshTWwUcJGb9U0RQt+s9Hn9pJy0R/ORftx8WKl+aUc2U0UM4a+Iwn6wvwlNbZtbvSWL6ycOYfvIwIoIrzhr/Wag8sHoHv9rYCBvz51euy51M9UnH85cXn+oT9cfIPy0z6xdah8qVZ09g9eZGXt/9IY+/vJuHa/NfbbR++z4uml5NZWWlp746yEFiZv2OJE4fN5TTxw0lIrjhC++zenMjz7y+h2dff4ff1O8FYHX9O5x36igHSjt8jsTMLNPc3MzS1dvY/cHHPPD8Dj44kP8azz85cwyzpozgvCkjmTbmpH7zYEmfIzEzO0aVlZXccMGU/NTXWRN4Ycs7rN32Lr946S1++cpuThwo/uyCKVxx1rh+FSjt8RGJmdlRHDpKqdv9AQ+vaaAF+k2g+PLfAg4SM0vV0tLCqld385u6Rh5as4MW4KRBFXzrsumcP2Ukp445sc8FioOkgIPEzEqlMFAeXruDloChx1Vwwxem8M3Zp/apk/I+R2Jm1gUqKiqYM3Mcs88Yw4XTR/PSjnf5wTOb+d7T9bzxzkfcfPHUPj3dVYyDxMysEyoqKrhkxlgunj6aP5owlF+s38Xyl3bx3Ou7WXDBFP7y4r51dHI0DhIzswSHjlAunj6ayaNO4MHV2/i7VfXs+eAAF00fzUXTR/f5R6/07dGZmXWTyspK/mrOaTz8tXO54bxJPLx2B3/xwDr+btUmWlpa2m+gF3OQmJmVyKE75m+/cia3XzmTEYMH8P1V9dy+fCPNzc3l7l6X8dSWmVmJVVRU8NXzTuHcySP459XbeGD1tvyjWM4/pU9eJuwgMTPrApI4bexJ3HbFDEA8tGYbT27YxTcvOY2vzJrUp86b9J2RmJn1QJWVldx+5QxuvzL/jY63LdvAQ2u205fu4UsKEkkjJK2UtCl7H95GvQVZnU2SFhSUnyPpZUn1kr6n7HivrXYlXSjpPUkvZq/vpPTfzKw7VFRUcP25p/Ddaz7HyCEDufeZep5+7e0+EyapRyS3AKsiYhqwKls/jKQRwG3AucAs4LaCwPkRsBCYlr3mdqDd30bE2dnrjsT+m5l1C0nMPmMM373mcwjxnZ+/0mfCJDVI5gH3Z8v3A1cVqXMZsDIi9kbEPmAlMFfSOOCkiHg+8j/JBwq270i7Zma9iiQunj6aO646E6DPhElqkIyJiF0A2fvoInXGAzsK1huysvHZcuvy9to9X9J6SU9ImtlWxyQtlFQrqbaxsfFYx2Vm1iX6Ypi0GySSnpL0SpHXvA7uo9h1bnGU8qP5HXBKRJwFfB/4eVsVI2JJROQiIlddXd3BrpqZdb3WYXLbsg3Uv/1hmXvVee0GSURcEhFnFnktA3ZnU1Rk728XaaIBmFiwPgHYmZVPKFJOW+1GxPsR8WG2/CtggKRRxzBeM7MeoTBMqipgzdY9vfYO+NSpreXAoauwFgDLitRZAVwqaXh2kv1SYEU2ZfWBpPOyq7VuKNi+aLuSxhZc2TUr6/+exDGYmZXFoTD52hen8j9+sZGH1+zolVNcqTck3gU8KulGYDtwLYCkHHBTRPxFROyVdCewNtvmjojYmy1/HfgpcDzwRPZqs13gy8DXJTUBHwPzozf+1M3MMpK4btZEIoJ//O1mxg49jounj+5Vd7/7i63MzHqA+t0fsOC+NUjiJ//xj5k25sRyd6nDX2zlO9vNzHqAqaOHcMdVZzKgUuzYu79XTXH5WVtmZj3AofMlAHc+vpGJIwb3iKOSjvARiZlZDyGJU0YMhqBXHZU4SMzMepCpo4fw36+Yyd1PvsYzdY29Ikw8tWVm1oNI4qLTq4kI7vzlBiYOz/X4KS4fkZiZ9TCHpriEij4CpKdxkJiZ9UBTRw/hH/7DOQA9fnrLQWJm1gMduiHxaw/W9vjncDlIzMx6KEGvmN5ykJiZ9VC9ZXrLQWJm1kNJQhLfeOh3PfpSYAeJmVkPNrX6BL59+Rn8zZOvsbnxo3J3pyjfR2Jm1oMduq9k0ojBTK0+odzdKcpBYmbWw0ni1NFDyt2NNnlqy8zMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLElSkEgaIWmlpE3Z+/A26i3I6myStKCg/BxJL0uql/Q9Zd8tKelaSRsktUjKtWrr1qx+naTLUvpvZmbpUo9IbgFWRcQ0YFW2fhhJI4DbgHOBWcBtBYHzI2AhMC17zc3KXwGuAZ5r1dYMYD4wM6t7r6TKxDGYmVmC1CCZB9yfLd8PXFWkzmXAyojYGxH7gJXAXEnjgJMi4vnIf+3XA4e2j4hXI6Kujf09EhEHImIrUE8+nMzMrExSg2RMROwCyN5HF6kzHthRsN6QlY3PlluXH01bbR1B0kJJtZJqGxsb22nWzMw6q90vtpL0FDC2yEeLOrgPFSmLo5R3pq0jCyOWAEsAcrlcz/yiYzOzPqDdIImIS9r6TNJuSeMiYlc2VfV2kWoNwIUF6xOAZ7PyCa3Kd7bTnQZg4jFuY2ZmXSh1ams5cOgqrAXAsiJ1VgCXShqenWS/FFiRTYV9IOm87GqtG9rYvvX+5ksaJGky+RP0axLHYGZmCVKD5C5gjqRNwJxsHUk5ST8GiIi9wJ3A2ux1R1YG8HXgx+RPmm8Gnsi2v1pSA3A+8LikFVlbG4BHgY3Ak8A3IqI5cQxmZpZA+Qum+rZcLhe1tbXl7oaZWa8iaV1E5Nqr5zvbzcwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySJAWJpBGSVkralL0Pb6PegqzOJkkLCsrPkfSypHpJ35OkrPxaSRsktUjKFdSvkfSxpBez19+n9N/MzNKlHpHcAqyKiGnAqmz9MJJGALcB5wKzgNsKAudHwEJgWvaam5W/AlwDPFdkn5sj4uzsdVNi/83MLFFqkMwD7s+W7weuKlLnMmBlROyNiH3ASmCupHHASRHxfEQE8MCh7SPi1YioS+ybmZl1g9QgGRMRuwCy99FF6owHdhSsN2Rl47Pl1uXtmSzp95J+I+nfdq7bZmZWKlXtVZD0FDC2yEeLOrgPFSmLo5QfzS5gUkTskXQO8HNJMyPi/SN2Ki0kP23GpEmTOthVMzM7Vu0GSURc0tZnknZLGhcRu7KpqreLVGsALixYnwA8m5VPaFW+s52+HAAOZMvrJG0GTgNqi9RdAiwByOVy7QWUmZl1UurU1nLg0FVYC4BlReqsAC6VNDw7yX4psCKbCvtA0nnZ1Vo3tLH9ZyRVS6rMlqeQP0G/JXEMZmaWIDVI7gLmSNoEzMnWkZST9GOAiNgL3AmszV53ZGUAXwd+DNQDm4Ensu2vltQAnA88LmlFVv+LwEuS1gOPATcVtGVmZmWg/AVTfVsul4va2iNmv8zM7CgkrYuIXHv1fGe7mZklcZCYmVkSB4mZmSVxkJiZWRIHiZmZJXGQmJlZEgeJmZklcZCYmVkSB4mZmSVxkJiZWRIHiZmZJXGQmJlZEgeJmZklcZCYmVkSB4mZmSVxkJiZWRIHiZmZJXGQmJlZEgeJmZklcZCYmVkSB4mZmSVxkJiZWRIHiZmZJXGQmJlZEkVEufvQ5SQ1Ats6ufko4J0Sdqc36G9j9nj7vv425lKN95SIqG6vUr8IkhSSaiMiV+5+dKf+NmaPt+/rb2Pu7vF6asvMzJI4SMzMLImDpH1Lyt2BMuhvY/Z4+77+NuZuHa/PkZiZWRIfkZiZWRIHSQdIul3Sm5JezF5fKnefuoOkb0kKSaPK3ZeuJulOSS9lf76/lnRyufvUlST9raTXsjH/i6Rh5e5TV5J0raQNklok9emrtyTNlVQnqV7SLd2xTwdJx90TEWdnr1+VuzNdTdJEYA6wvdx96SZ/GxGfi4izgV8C3yl3h7rYSuDMiPgc8Dpwa5n709VeAa4Bnit3R7qSpErgh8DlwAzgOkkzunq/DhJryz3AfwX6xUm0iHi/YPUE+vi4I+LXEdGUra4GJpSzP10tIl6NiLpy96MbzALqI2JLRBwEHgHmdfVOHSQdd3M2DXCfpOHl7kxXknQl8GZErC93X7qTpMWSdgDX0/ePSAr9OfBEuTthJTEe2FGw3pCVdamqrt5BbyHpKWBskY8WAT8C7iT/W+qdwP8i/4+v12pnvP8NuLR7e9T1jjbmiFgWEYuARZJuBW4GbuvWDpZYe+PN6iwCmoCl3dm3rtCR8fYDKlLW5UfXDpJMRFzSkXqS/pH8HHqv1tZ4Jf0RMBlYLwnyUx6/kzQrIt7qxi6WXEf/jIGHgMfp5UHS3nglLQD+BJgdfeA+gGP48+3LGoCJBesTgJ1dvVNPbXWApHEFq1eTP3HXJ0XEyxExOiJqIqKG/F/Mz/f2EGmPpGkFq1cCr5WrL91B0lzg28CVEbG/3P2xklkLTJM0WdJAYD6wvKt36iOSjvkbSWeTP0R8A/hP5e2OdYG7JJ0OtJB/UvRNZe5PV/sBMAhYmR15ro6IPjtmSVcD3weqgcclvRgRl5W5WyUXEU2SbgZWAJXAfRGxoav36zvbzcwsiae2zMwsiYPEzMySOEjMzCyJg8TMzJI4SMzM+pjsCRxvS2r3VgVJp0halT2541lJx/y4HAeJmVnf81Ngbgfr/k/ggewBnncA3z3WnTlIzMz6mIh4DthbWCZpqqQnJa2T9FtJ07OPZgCrsuVn6MRDHh0kZmb9wxLgP0fEOcC3gHuz8vXAn2bLVwMnShp5LA37znYzsz5O0hDgC8D/zZ5kAPknG0A+VH4g6c/If1/Lm+Qf5NlhDhIzs76vAng3++K2w0TETvJf+nUocP40It471sbNzKwPy764baukawGUd1a2PErSoSy4FbjvWNt3kJiZ9TGSHgaeB06X1CDpRvJf2HajpPXABv5wUv1CoE7S68AYYPEx788PbTQzsxQ+IjEzsyQOEjMzS+IgMTOzJA4SMzNL4iAxM7MkDhIzM0viIDEzsyQOEjMzS/L/AXA5a2m7U0JDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.scatter(x,y, s=0.2)\n",
    "plt.scatter(0,0, c='black')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### It can be seen that as the particle approaches the massive body, it acquires axial velocity due to spin and frame-dragging effect of the body."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another Way of calculating Geodesics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 0.3 * u.m\n",
    "Attractor = Body(name=\"BH\", mass=1.989e30 * u.kg, a=a)\n",
    "init_conditions = BoyerLindquistDifferential(49.95e5 * u.km, np.pi / 2 * u.rad, \n",
    "                                             np.pi * u.rad, 0 * u.km / u.s, \n",
    "                                             0 * u.rad / u.s, 0 * u.rad / u.s,\n",
    "                                             a)\n",
    "Particle = Body(differential=init_conditions, parent=Attractor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "geodesic = Geodesic(body=Particle, time=0 * u.s, end_lambda=33932.90, \n",
    "                    step_size=1.2, metric=Kerr)\n",
    "ans = geodesic.trajectory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = ans[:,1], ans[:,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEJCAYAAAC+I6F6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGlRJREFUeJzt3X2Q1dWd5/H3p7sBBeS5eZAHGxBFcKIb76DG2pSKKGZHUSfWYszK7DhhTcWdTG1lNrrWBlfLis7srjt5MDNMxkQd1LhWTSAxShA1pnZEaBJRAVsaEGhBbAEfUbC7v/vH/WEuzW266XO7b3fzeVXdur/fued3fuc0D5/+nd/DVURgZmbWWRXl7oCZmfVuDhIzM0viIDEzsyQOEjMzS+IgMTOzJA4SMzNL4iAxM7MkDhIzM0viIDEzsyRV5e5Adxg1alTU1NSUuxtmZr3K2rVr34mI6vbqHRdBUlNTQ21tbbm7YWbWq0ja1pF6ntoyM7MkDhIzM0viIDEzsyQOEjMzS+IgMTOzJA4SMzNL4iAxM7MkJQkSSXMl1Umql3RLkc8HSPpZ9vmLkmoKPrs1K6+TdFlB+f2S3pb0aqu2RkhaIWlT9j68FGMwM7POSQ4SSZXAD4HLgRnAdZJmtKp2I7AvIk4F7gXuybadAcwHZgJzgfuy9gB+mpW1dguwMiKmASuzdTMzK5NSHJHMAuojYktEHAQeBea1qjMPeCBbfhyYLUlZ+aMRcSAitgL1WXtExPPA3iL7K2zrAeCqEozBzMw6qRRBMh7YUbDekJUVrRMRTcB7wMgObtvamIjYlbW1CxhdrJKkhZJqJdU2NjZ2cChmZnasShEkKlIWHazTkW07JSIWR0QuInLV1e0+c8zMzDqpFEHSAEwsWJ8A7GyrjqQqYCj5aauObNvabknjsrbGAW93uudmZpasFEGyBpgmabKk/uRPni9rVWcZsCBb/jLwTEREVj4/u6prMjANWN3O/grbWgAsLcEYzMysk5KDJDvncTOwHNgIPBYR6yXdIenKrNo/ASMl1QP/hexKq4hYDzwGbACeAr4REc0Akh4BXgBOl9Qg6casrbuBOZI2AXOydTMzKxPlDwz6tlwuF/4+EjOzYyNpbUTk2qvnO9vNzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJI4SMzMLImDxMzMkjhIzMwsiYPEzMySOEjMzCyJg8TMzJKUJEgkzZVUJ6le0i1FPh8g6WfZ5y9Kqin47NasvE7SZe21KemnkrZKeil7nV2KMZiZWedUpTYgqRL4ITAHaADWSFoWERsKqt0I7IuIUyXNB+4B/r2kGcB8YCZwMvC0pNOybY7W5l9HxOOpfTczs3SlOCKZBdRHxJaIOAg8CsxrVWce8EC2/DgwW5Ky8kcj4kBEbAXqs/Y60qaZmfUApQiS8cCOgvWGrKxonYhoAt4DRh5l2/bavEvSy5LulTSgBGMwM7NOKkWQqEhZdLDOsZYD3ApMB/4YGAF8u2inpIWSaiXVNjY2FqtiZmYlUIogaQAmFqxPAHa2VUdSFTAU2HuUbdtsMyJ2Rd4B4Cfkp8GOEBGLIyIXEbnq6upODs3MzNpTiiBZA0yTNFlSf/Inz5e1qrMMWJAtfxl4JiIiK5+fXdU1GZgGrD5am5LGZe8CrgJeLcEYzMysk5Kv2oqIJkk3A8uBSuD+iFgv6Q6gNiKWAf8EPCSpnvyRyPxs2/WSHgM2AE3ANyKiGaBYm9kul0iqJj/99RJwU+oYzMys85Q/MOjbcrlc1NbWlrsbZma9iqS1EZFrr57vbDczsyQOEjMzS+IgMTOzJA4SMzNL4iAxM7MkDhIzM0viIDEzsyQOEjMzS+IgMTOzJA4SMzNL4iAxM7MkDhIzM0viIDEzsyQOEjMzS+IgMTOzJA4SMzNL4iAxM7MkDhIzM0viIDEzsyQOEjMzS+IgMTOzJA4SMzNL4iAxM7MkDhIzM0viIDEzsyQOEjMzS1KSIJE0V1KdpHpJtxT5fICkn2WfvyippuCzW7PyOkmXtdempMlZG5uyNvuXYgxmZtY5yUEiqRL4IXA5MAO4TtKMVtVuBPZFxKnAvcA92bYzgPnATGAucJ+kynbavAe4NyKmAfuyts3MbMkSqKmBior8+5Il3bLbUhyRzALqI2JLRBwEHgXmtaozD3ggW34cmC1JWfmjEXEgIrYC9Vl7RdvMtrk4a4OszatKMAYzs95tyRJYuBC2bYOI/PvChd0SJlUlaGM8sKNgvQE4t606EdEk6T1gZFa+qtW247PlYm2OBN6NiKYi9btEzS1PdGXzPV4F0AJUFpQpW28ueK8q8nlTtlwlkPJ/t5sC+lVAv0r4tDmrL6iqgKaWfPnRSH/Ytn9Vvs1Dy23V718lDnwany0fbAoG9Mv/DhURh63/YTtxQlUFH3/afMTyif0qiQg+aWrhxKzDrdclMbBfBR8dbAHisLL9nwYD++mw90H9K4iA/Z+2MKh/RfaTO9rPQQzqX8mHB5o+W/7oYPMR74MHVBERny0fTWVlJedNrear59dQWdnOH4T1PLfdBvv3H162f3++/Prru3TXpQiSYn/jo4N12iovdqR0tPpHdkpaCCwEmDRpUrEq1gEt2Xtzq/KmNt6LCg77U/q4paDhQ5pbvXfEgTaWi3ag9XLrDrReb92ZjiwXW+9dnty4h4dffIOzJwxl9NBBfPOSaVRVleK/Cety27cfW3kJleJvSAMwsWB9ArCzjToNkqqAocDedrYtVv4OMExSVXZUUmxfAETEYmAxQC6XKxo21j4fkRw/RyRvvneQ9bv3U9f4CXWNnwDwi5d2MKtmOP369ef2K2fQr1+/o7ZhZTRpUn46q1h5FytFkKwBpkmaDLxJ/uT5V1rVWQYsAF4Avgw8ExEhaRnwsKT/DZwMTANWk/9XdESb2TbPZm08mrW5tARjaNMbd/+7rmzerMdobm7mwf+3hRe37uGjA838dsu7vPHup7zx0tsAPLvxLS6YOpwPDwbnTR3lKbCe5q678udECqe3Bg7Ml3cxRaT/si7pS8D/If+L6P0RcZekO4DaiFgm6QTgIeDfkD8SmR8RW7JtbwP+nPwvsH8VEU+21WZWPoV8iIwAfg98NSKOOrGRy+WitrY2eZxmx4uWlhZ+/epOlv7+TSJaeLpuL02t/qs4a9xAxg45gWED+3karKdYsiR/TmT79vyRyF13JZ0fkbQ2InLt1itFkPR0DhKzzosINr65jx88s4lB/cSWPZ+wtuGDI+rVDOvPrJphgKisrOTiM8Ywe8ZYKip833Nv5SAp4CAxK52WlhZWrN/FMxve4t39B6moEE/X7eHTItcrXFAzlPHDTvhsvaqqihu+UMPp44aSv5rfejIHSQEHiVnXKTxiGZxdKLB178es2f5+0foCLps+kpqRA9m7/yBCVA85kSvOPpnTxjpgehIHSQEHiVn3am5u5qF/3crGXe8RLX/4P2bD7o94ddeHRbepBOZMH8mU6sGMHT6Q6889xSfzy8xBUsBBYtYzNDU1ce+v63jnw08YfmK/z45Ituz9mDXb3jus7gU1Q5kw4kRGDhrgYCkTB0kBB4lZz1YYMB8eaOaJ9Y1H1Pn8hJM4t2Y4KH+u5ZuzfZVYV3OQFHCQmPUeLS0tPL3hLV7atpc9Hx2gQhX8y7pdfNLq+uPPjx/MtDGDGTX4BM4+ZTizz/AVYqXmICngIDHrvSKC13a+y0MvvMHQE6pAULv9fdZse/ewepfPqKZm5Ik+WikhB0kBB4lZ39LeVNgXpw7nzPFDABg7zOdXOqujQeLINrNep6qqir/+0kwgPxU2b8NbrNu+j5ZoYeNbH/Hcpj08v3nfZ/WfXr+bi6aPoqqqiq/MmuRQKTEfkZhZn9Lc3Mw/v/AGu9//mH37P2XpurfYX3C35BenjnCodJCntgo4SMyOTxFB3a73eOLlnbRE8JN/3X5YqHxpRjVTRg/mrInDfLK+CE9tmdlxTxLTTx7G9JOHERFccdb4z0LlwVU7+NWGRtiQP79yXe5kqoecyF9efKpP1B8j/7TM7LjQOlSuPHsCqzY38vruD3nild08Upv/aqN12/dx0fRqKisrPfXVQQ4SMzvuSOL0cUM5fdxQIoIbvvA+qzY38uzre3ju9Xf4Tf1eAFbVv8N5p45yoLTD50jMzDLNzc0sWbWN3R98zIMv7OCDA/mv8fyTM8cwa8oIzpsykmljhhw3D5b0ORIzs2NUWVnJDRdMyU99nTWBF7e8w5pt7/KLl9/il6/u5qT+4s8umMIVZ407rgKlPT4iMTM7ikNHKXW7P+CR1Q20wHETKL78t4CDxMxStbS0sHLjbn5T18jDq3fQAgwZUMG3LpvO+VNGcuqYk/pcoDhICjhIzKxUCgPlkTU7aAkYekIFN3xhCt+cfWqfOinvcyRmZl2goqKCOTPHMfuMMVw4fTQv73iXHzy7me89U88b73zEzRdP7dPTXcU4SMzMOqGiooJLZozl4umj+aMJQ/nFul0se3kXz7++mwUXTOEvL+5bRydH4yAxM0tw6Ajl4umjmTxqEA+t2sbfraxnzwcHuGj6aC6aPrrPP3qlb4/OzKybVFZW8ldzTuORr53LDedN4pE1O/iLB9fydys30dLS0n4DvZiDxMysRA7dMX/7lTO5/cqZjBjYj++vrOf2ZRtobm4ud/e6jKe2zMxKrKKigq+edwrnTh7BP6/axoOrtuUfxXL+KX3yMmEHiZlZF5DEaWOHsOiKGYB4ePU2nlq/i29echpfmTWpT5036TsjMTPrgSorK7n9yhncfmX+Gx0XLV3Pw6u305fu4UsKEkkjJK2QtCl7H95GvQVZnU2SFhSUnyPpFUn1kr6n7HivrXYlXSjpPUkvZa/vpPTfzKw7VFRUcP25p/Ddaz7HyMH9ue/Zep557e0+EyapRyS3ACsjYhqwMls/jKQRwCLgXGAWsKggcH4ELASmZa+5HWj3txFxdva6I7H/ZmbdQhKzzxjDd6/5HEJ85+ev9pkwSQ2SecAD2fIDwFVF6lwGrIiIvRGxD1gBzJU0DhgSES9E/if5YMH2HWnXzKxXkcTF00dzx1VnAvSZMEkNkjERsQsgex9dpM54YEfBekNWNj5bbl3eXrvnS1on6UlJM9vqmKSFkmol1TY2Nh7ruMzMukRfDJN2g0TS05JeLfKa18F9FLvOLY5SfjS/A06JiLOA7wM/b6tiRCyOiFxE5KqrqzvYVTOzrtc6TBYtXU/92x+WuVed126QRMQlEXFmkddSYHc2RUX2/naRJhqAiQXrE4CdWfmEIuW01W5EvB8RH2bLvwL6SRp1DOM1M+sRCsOkqgJWb93Ta++AT53aWgYcugprAbC0SJ3lwKWShmcn2S8FlmdTVh9IOi+7WuuGgu2LtitpbMGVXbOy/u9JHIOZWVkcCpOvfXEq/+MXG3hk9Y5eOcWVekPi3cBjkm4EtgPXAkjKATdFxF9ExF5JdwJrsm3uiIi92fLXgZ8CJwJPZq822wW+DHxdUhPwMTA/euNP3cwsI4nrZk0kIvjH325m7NATuHj66F5197u/2MrMrAeo3/0BC+5fjSR+8h//mGljTip3lzr8xVa+s93MrAeYOnowd1x1Jv0qxY69+3vVFJeftWVm1gMcOl8CcOcTG5g4YmCPOCrpCB+RmJn1EJI4ZcRACHrVUYmDxMysB5k6ejD//YqZ3PPUazxb19grwsRTW2ZmPYgkLjq9mojgzl+uZ+LwXI+f4vIRiZlZD3Noikuo6CNAehoHiZlZDzR19GD+4T+cA9Djp7ccJGZmPdChGxK/9lBtj38Ol4PEzKyHEvSK6S0HiZlZD9VbprccJGZmPZQkJPGNh3/Xoy8FdpCYmfVgU6sH8e3Lz+BvnnqNzY0flbs7Rfk+EjOzHuzQfSWTRgxkavWgcnenKAeJmVkPJ4lTRw8udzfa5KktMzNL4iAxM7MkDhIzM0viIDEzsyQOEjMzS+IgMTOzJA4SMzNL4iAxM7MkDhIzM0viIDEzsyQOEjMzS+IgMTOzJElBImmEpBWSNmXvw9uotyCrs0nSgoLycyS9Iqle0veUfbekpGslrZfUIinXqq1bs/p1ki5L6b+ZmaVLPSK5BVgZEdOAldn6YSSNABYB5wKzgEUFgfMjYCEwLXvNzcpfBa4Bnm/V1gxgPjAzq3ufpMrEMZiZWYLUIJkHPJAtPwBcVaTOZcCKiNgbEfuAFcBcSeOAIRHxQuS/9uvBQ9tHxMaIqGtjf49GxIGI2ArUkw8nMzMrk9QgGRMRuwCy99FF6owHdhSsN2Rl47Pl1uVH01ZbR5C0UFKtpNrGxsZ2mjUzs85q94utJD0NjC3y0W0d3IeKlMVRyjvT1pGFEYuBxQC5XK5nftGxmVkf0G6QRMQlbX0mabekcRGxK5uqertItQbgwoL1CcBzWfmEVuU72+lOAzDxGLcxM7MulDq1tQw4dBXWAmBpkTrLgUslDc9Osl8KLM+mwj6QdF52tdYNbWzfen/zJQ2QNJn8CfrViWMwM7MEqUFyNzBH0iZgTraOpJykHwNExF7gTmBN9rojKwP4OvBj8ifNNwNPZttfLakBOB94QtLyrK31wGPABuAp4BsR0Zw4BjMzS6D8BVN9Wy6Xi9ra2nJ3w8ysV5G0NiJy7dXzne1mZpbEQWJmZkkcJGZmlsRBYmZmSRwkZmaWxEFiZmZJHCRmZpbEQWJmZkkcJGZmlsRBYmZmSRwkZmaWxEFiZmZJHCRmZpbEQWJmZkkcJGZmlsRBYmZmSRwkZmaWxEFiZmZJHCRmZpbEQWJmZkkcJGZmlsRBYmZmSRwkZmaWxEFiZmZJHCRmZpbEQWJmZkmSgkTSCEkrJG3K3oe3UW9BVmeTpAUF5edIekVSvaTvSVJWfq2k9ZJaJOUK6tdI+ljSS9nr71P6b2Zm6VKPSG4BVkbENGBltn4YSSOARcC5wCxgUUHg/AhYCEzLXnOz8leBa4Dni+xzc0Scnb1uSuy/mZklSg2SecAD2fIDwFVF6lwGrIiIvRGxD1gBzJU0DhgSES9ERAAPHto+IjZGRF1i38zMrBukBsmYiNgFkL2PLlJnPLCjYL0hKxufLbcub89kSb+X9BtJ/7Zz3TYzs1Kpaq+CpKeBsUU+uq2D+1CRsjhK+dHsAiZFxB5J5wA/lzQzIt4/YqfSQvLTZkyaNKmDXTUzs2PVbpBExCVtfSZpt6RxEbErm6p6u0i1BuDCgvUJwHNZ+YRW5Tvb6csB4EC2vFbSZuA0oLZI3cXAYoBcLtdeQJmZWSelTm0tAw5dhbUAWFqkznLgUknDs5PslwLLs6mwDySdl12tdUMb239GUrWkymx5CvkT9FsSx2BmZglSg+RuYI6kTcCcbB1JOUk/BoiIvcCdwJrsdUdWBvB14MdAPbAZeDLb/mpJDcD5wBOSlmf1vwi8LGkd8DhwU0FbZmZWBspfMNW35XK5qK09YvbLzMyOQtLaiMi1V893tpuZWRIHiZmZJXGQmJlZEgeJmZklcZCYmVkSB4mZmSVxkJiZWRIHiZmZJXGQmJlZEgeJmZklcZCYmVkSB4mZmSVxkJiZWRIHiZmZJXGQmJlZEgeJmZklcZCYmVkSB4mZmSVxkJiZWRIHiZmZJXGQmJlZEgeJmZklcZCYmVkSB4mZmSVRRJS7D11OUiOwrZObjwLeKWF3eoPjbcweb993vI25VOM9JSKq26t0XARJCkm1EZErdz+60/E2Zo+37zvextzd4/XUlpmZJXGQmJlZEgdJ+xaXuwNlcLyN2ePt+463MXfreH2OxMzMkviIxMzMkjhIOkDS7ZLelPRS9vpSufvUHSR9S1JIGlXuvnQ1SXdKejn78/21pJPL3aeuJOlvJb2WjflfJA0rd5+6kqRrJa2X1CKpT1+9JWmupDpJ9ZJu6Y59Okg67t6IODt7/arcnelqkiYCc4Dt5e5LN/nbiPhcRJwN/BL4Trk71MVWAGdGxOeA14Fby9yfrvYqcA3wfLk70pUkVQI/BC4HZgDXSZrR1ft1kFhb7gX+K3BcnESLiPcLVgfRx8cdEb+OiKZsdRUwoZz96WoRsTEi6srdj24wC6iPiC0RcRB4FJjX1Tt1kHTczdk0wP2Shpe7M11J0pXAmxGxrtx96U6S7pK0A7ievn9EUujPgSfL3QkrifHAjoL1hqysS1V19Q56C0lPA2OLfHQb8CPgTvK/pd4J/C/y//h6rXbG+9+AS7u3R13vaGOOiKURcRtwm6RbgZuBRd3awRJrb7xZnduAJmBJd/atK3RkvMcBFSnr8qNrB0kmIi7pSD1J/0h+Dr1Xa2u8kv4ImAyskwT5KY/fSZoVEW91YxdLrqN/xsDDwBP08iBpb7ySFgB/AsyOPnAfwDH8+fZlDcDEgvUJwM6u3qmntjpA0riC1avJn7jrkyLilYgYHRE1EVFD/i/m53t7iLRH0rSC1SuB18rVl+4gaS7wbeDKiNhf7v5YyawBpkmaLKk/MB9Y1tU79RFJx/yNpLPJHyK+Afyn8nbHusDdkk4HWsg/KfqmMvenq/0AGACsyI48V0VEnx2zpKuB7wPVwBOSXoqIy8rcrZKLiCZJNwPLgUrg/ohY39X79Z3tZmaWxFNbZmaWxEFiZmZJHCRmZpbEQWJmZkkcJGZmfUz2BI63JbV7q4KkUyStzJ7c8ZykY35cjoPEzKzv+Skwt4N1/yfwYPYAzzuA7x7rzhwkZmZ9TEQ8D+wtLJM0VdJTktZK+q2k6dlHM4CV2fKzdOIhjw4SM7Pjw2LgP0fEOcC3gPuy8nXAn2bLVwMnSRp5LA37znYzsz5O0mDgC8D/zZ5kAPknG0A+VH4g6c/If1/Lm+Qf5NlhDhIzs76vAng3++K2w0TETvJf+nUocP40It471sbNzKwPy764baukawGUd1a2PErSoSy4Fbj/WNt3kJiZ9TGSHgFeAE6X1CDpRvJf2HajpHXAev5wUv1CoE7S68AY4K5j3p8f2mhmZil8RGJmZkkcJGZmlsRBYmZmSRwkZmaWxEFiZmZJHCRmZpbEQWJmZkkcJGZmluT/A0DTWXhqD1XLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "plt.scatter(x,y, s=0.2)\n",
    "plt.scatter(0,0, c='red')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}